#Alexander F. Gazmararian
#afg2@princeton.edu
#January 19, 2024

#Load packages
library(tidyverse)
library(sandwich)
library(lmtest)
library(modelsummary)
library(kableExtra)
library(broom)
library(scales)
library(viridis)
library(here)

#Load functions
source(here("code", "fun", "savefig.r"))
source(here("code", "fun", "gofmap.r"))
source(here("code", "fun", "fix_txt.r"))

#Load data
g <- readRDS(here("data", "qualtrics", "natsurvey.rds"))

#Load baseline model specification
f.base <- y ~ age + female + black + hispanic + aapi + rural_bin + college + fullemploy + pid3 + ownhome + income5 + buyself_bin + suitable_bin + index_trust + altruism + risk

#Create coefficient mapping for tables
coefmap <- c(
  "poweroffer"="Utility Offers Program",
  "poweroffer:risk"="Utility Offers Program x Risk Preferences",
  "poweroffer:altruism"="Utility Offers Program x Altruism",
  "poweroffer:index_trust" = "Utility Offers Program x Trust Index",
  "poweroffer:buyself_bin"="Utility Offers Program x Solar Interest",
  "poweroffer:pid3Republican"="Utility Offers Program x Republican",
  "poweroffer:pid3Neither"="Utility Offers Program x Neither"
)

#Estimate models
m <- list()
# Donation allocation outcome
m[[1]] <- lm(update(f.base, power_num ~ poweroffer + .), g)
m[[2]] <- lm(update(f.base, power_num ~ poweroffer * risk + .), g)
m[[3]] <- lm(update(f.base, power_num ~ poweroffer * altruism + .), g)
m[[4]] <- lm(update(f.base, power_num ~ poweroffer * index_trust + .), g)
m[[5]] <- lm(update(f.base, power_num ~ poweroffer * buyself_bin + .), g)
# Cares about the community outcome
m[[6]] <- lm(update(f.base, powercom_num ~ poweroffer + .), g)
m[[7]] <- lm(update(f.base, powercom_num ~ poweroffer * risk + .), g)
m[[8]] <- lm(update(f.base, powercom_num ~ poweroffer * altruism + .), g)
m[[9]] <- lm(update(f.base, powercom_num ~ poweroffer * index_trust + .), g)
m[[10]] <- lm(update(f.base, powercom_num ~ poweroffer * buyself_bin + .), g)
# Create table
file <- here("output", "tables", "tab_power.txt")
modelsummary(
  m,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC2",
  coef_map = coefmap,
  gof_map = gof_map,
  gof_omit = "R2|IC|Log|F|RMSE|Err",
  escape = FALSE,
  output = "latex",
  add_rows = data.frame(
    c("Covariates"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes"),
    c("Yes")
  )
) %>%
  add_header_above(c(" " = 1, "Favorability" = 5, "Cares About Community" = 5)) %>%
  kable_styling(latex_options=c("HOLD_position")) %>%
  cat(., file = file)
fix_txt(file)
